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Abstract 

We consider stochastic volatility models using piecewise constant pa- 
rameters. We suggest a hybrid optimization algorithm for fitting the mod- 
els to a volatility surface and provide some numerical results. Finally, we 
provide an outlook on how to further improve the calibration procedure. 



1 Introduction 

Over the past two decades various stochastic volatility models have been intro- 
duced to explain stylized facts observed in the market. Some of the most popular 
models are the Heston stochastic volatility model and an extension allowing for 
jumps, cf. Heston [1993 and Bates 1996 



model introduced by Madan, Carr, and Chang 1998 



respectively, the Variance-Gamma 
the model introduced by 



sec 



Carr, Geman, Madan, and Yor 2003 



Barndorff-Nielsen and Shephard 2001 and Levy models with stochastic time, 



There are basically two approaches to calibrate a stochastic volatility model. 
It is possible to calibrate the model to historical data or to (liquid) present op- 
tion prices. For the first approach there exist various methods like maximum 
likelihood methods, efficient method of moments, or Markov Chain Monte Carlo 
methods, see A'it-Sahalia and Kimmel 2007 , Gallant and Tauchen 1996 and 



Friihwirth-Schnatter and Soegner 2007 , respectively, and the references therein. 



For various applications the second approach is favorable since the fitted model 
explains current prices observed on the market, hence we will subsequently con- 
sider the latter approach. 

For fitting the model we use a hybrid optimization algorithm consisting of 
a genetic algorithm to roughly locate the minimum and subsequently a pattern 
search algorithm with the objective to refine the result. This algorithm allows 
to take nonlinear constraints like the Feller condition for the Heston model into 
account. Evolutionary algorithms like the genetic algorithm have already been 



used for calibration e.g. by Hamida and Cont 2005 



*With best thanks to Cedric de Masson d'Autume who contributed some of the code used 
in the numerical example. 

twolfgang.putschoegl@unicreditgroup.at, Counterparty Risk Management & Analysis, 
UniCredit Bank Austria, A-1090 Vienna. 

Disclaimer: The opinion expressed here is that of the author and does not reflect the views 
or policies of his employer. 
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The outline of this paper is as follows. We start in Section [2] with an in- 
troduction to different models for option pricing. We outline in Section [3] how 
to fit a stochastic volatility model to a given volatility surface and give our 
conclusions in Section [H 



2 Model Setup 



Let X = {X t ) te i 0iT ] denote the stock price process and 4>t(u) the characteristic 
function of the log-price ln(JQ), defined by 

<f>t(u) = E ex.p(tu\n(X t ))\X 0l v . 

For ease of notation we assume the risk-free interest rate r t = r to be constant. 



2.1 The Heston Model 



Heston 1993 introduced a stochastic volatility model which satisfies a lot of 



stylized facts like fat tails, volatility clustering, skewness etc., cf. Cont 2001 



The Heston model is defined under the risk neutral measure by the coupled 
two-dimensional SDE 



dX t = rX t dt + yfoXt dW t , 
dv t — k(6 — v t ) dt + a^/v^dWt 



(2.1) 
(2.2) 



where W = (Wt)te[o,T\ an d W — {Wt)te[o,T\ are one-dimensional Brownian 
motions with correlation p and v = {vt)te[o,T] denotes the variance process. The 
involved parameters are the rate of mean reversion k > 0, the long term variance 
> 0, the volatility of variance a > (often called the vol of vol), the correlation 
— 1 < p < 1, and the initial variance Vq. If the Feller condition 2k9 > a 2 is 
satisfied, then the variance process is always positive and cannot reach zero. 
In absence of the stochastic factor, we have an exponential attraction to long 
term variance, the equilibrium point being v t = 9. Typically, the correlation p 
is negative, hence a down-move in the stock price is positively correlated with 
an up-move in the volatility. 

2.1.1 Time-Independent Parameters 

Note that there are two representations of the characteristic function. The first 
one used e.g. by Heston [1993] or |Kahl and Jackel 2005 reads 

<f>t(u) = exp(jw(lnX + rt)) 

x exp(6>Ker~ 2 ((K - paiu + d)t - 21n((l - gie dt )/(l - g{)))) 
x exp(v cr~ 2 (K - poiu + d)(l - e dt )/(l - gie dt j) , 



and the second one e.g. used by Schoutens, Simons, and Tistaert 2004 and 
Gatheral 2006] is given by 



4 2 \u) = exp(m{\nX + rt)) (2.3) 
x exp(6>Ker~ 2 ((K - paiu - d)t - 21n((l - g 2 e dt )/(l - g 2 )))) 
x exp(w cr _2 (K - paiu - d)(l - e dt )/(l - g2d dt )) , 
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where 



d 
9i 



[paui — k)^ 
k — paiu + d 



K - 



paiu 



a 2 (iu ■ 
92 -- 



K — paiu — d 1 
k — poiu + d g\ 



(2.4) 
(2.5) 



Albrecher, Mayer, Schoutens, and Tistaert 2007 prove in a detailed analysis 
that <fi^ L ) is unstable under certain conditions whereas (f>^ is stable under the 
full parameter space. 

2.1.2 Piecewise Constant Parameters 

To obtain a good fit for a volatility surface, i.e., for several maturities at once, is 
fairly difficult to obtain using a stationary process. One way to obtain a good fit 
is to use processes for the price process having piecewise constant parameters. 
The solution for the characteristic function for piecewise constant parameters is 



given by Mikhailov and Nogel [20031 who derived it using a Computer Algebra 



system. The full derivation of the solution is provided by Elices 2007 who also 



presents a formula which is in line with the findings of Albrecher et al. 
The solution is given by 



2007 



4>t(u) = exp(C + Dv + iuln(X )) , 



where 



C = 
D = 



+ ~ (~21n( 1 i gC ~" ) + (« - paiu - d)t) + C° , 

k — patu — d — D°a 2 



itru 

K — paui ( g 



ge~ dt 

-dt 



1 — ge 



paiu 



d - D a a 2 ' 



with d as in (2.4 1 and g = 52 as in (2.5) 



The solution is close to the Heston one ( |2.3[ ). The time interval to matu- 
rity [0, T] is divided into n sub-intervals [0, ti], . ■ ■ , [tfc, tfc+i], . . . , [t n -i, T] where 
i/c, k — 1, . . . , n — 1 is the time of model parameter jumps. Model parameters 
are constant during [t^, tk+i] but different for each sub interval. The initial con- 
dition for the first sub interval from the end [0, ti] is zero. For the second sub 
interval [£1, £2] we use C an d D from the first sub interval for C° and D°, the ini- 
tial conditions. The same procedure is repeated at each time £j., k = 2, . . . , n— 1 
of the parameters jumps. 

2.2 The Heston Model with Jumps 



Bates 1996 extended the Heston model by introducing jumps in the asset price. 



The coupled two-dimensional SDE (2.1)-(2.2) becomes 



dX t = (r- Xpj)X t dt + ^T t X t dW t + J t dN t , 
dv t = k{6 - v t ) dt + ay/v~ t dW t , 



with the same notations as in (2.1) and (2.2) and with N — (Ar t ) tg [ ,r] being 
an independent Poisson process with intensity A > (i.e., E[N t ] = Xt) and 
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J = (Jt)te[o,T] is the percentage jump size that is log-normally i.i.d. over time 
with mean fij. The standard deviation of ln(l + Jt) is <jj and 

ln(l + J t )^M(Hl + ^j)~^-,a 2 j) . 



2.2.1 Time-Independent Parameters 

Since the jumps and the diffusion part are independent, the characteristic func- 
tion for ln(Xj) is defined by 

cfi t (u) = <fi?(u)cfi J t {u) , 

where 4>? = 4>t as m (2.3 1 is the characteristic function for the diffusion part 
and <fi({u) the characteristic function corresponding to the jump part is given 

by 

cfi({u) = exp(-XfXjmt + Xt((l + /j,j) m exp(a 2 j(iu/2)(iu - 1)) -if) . 



2.2.2 Piecewise Constant Parameters 

We still define the characteris tic fun ction as 4>t{u) = 4>f ( u )4>t ( u ) with <fi^ cor- 
responding to (fit as in Section 2.1.2 and <fi~[ given by 



= E 



N-1 r t H +l 



exp 



fe=0 ' ,lk 



(fit (it) = E exp (iu j J s dN s 

N-1 

= [] E[exp(z 

N-1 

= [] e X p(-\W$\uA tk + A«A tfc ((l + M f )« 



J, dN f 



k=0 



x exp{af )2 {iu/2){iu- 1)) - 1)) , 



where and are the parameters of the jumps during the sub interval 
[tk,tk+i] and A ifc = t^+i — 4*. 



2.3 The Variance Gamma Process 

The Variance-Gamma (VG) process is obtained by evaluating Brownian motion 
(with constant drift 9 and volatility a) at a random time change given by a 
Gamma process. Let 

b t {6,a) = et + aW t , 

where W — (W / t ) t6 [o,T] is a Brownian motion. Then b t (9,a) is a Brownian 
motion with drift 9 and volatility a. 

The Gamma process 7t(/x, v) with mean rate /j. and variance rate v is a 
process of independent Gamma increments over time intervals [4,4 + h]. The 
VG process G t (cr, v, 9) is then defined by 

G t (a,v,0) = b lt{1 ^ } (9,a) . 
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The involved parameters are the volatility of the Brownian motion a, the vari- 
ance rate time change v and the drift in the Brownian motion 9. We can control 
the skewness of the distribution via 9 and the kurtosis with v. The resulting 
risk-neutral process for the stock price is 

X t = X a exp(rt + G t (a, 9, v) + ut) , 

where to = {1/v) ln(l — 9v — a 2 v/2). 

2.3.1 Time-Independent Parameters 

The characteristic function of \n(X t ) reads 

4> t (u) = exp(ln(X + (r + u)t)) (l - i9vu + a 2 u 2 v /2)~ t/v . 
From this formula, one can obtain a closed formula for both the density function 



and the option price, for details we refer to Madan et al. 1998 



2.3.2 Piecewise Constant Parameters 

For piecewise constant parameters we can use the same method as in Sec- 
tion 



2.2.2 For readability we use the notation G t instead of G t (<7, 9, u), AG tfc 



Gt k+l — Gt k and t/v = t. Note that Go = and that the increments AGt fc are 
independent. Then 4>t(u) is given by 

4>t{u) = E exp(iuln(X t )) 
= E exp^u(ln(X ) + (r + w)t)J exp(zuG t ) 

JV-1 

= exp(m(ln(Xo) + (r + to 

)*)) E [I1 exp(m(AG : 



fc=0 
N-l 



exp 



(iu(ln(X ) + {r + w)t)) Yl E[exp^u(AG t 

k=0 
N—l 

cxp^iu(ln(Xo) + (r + lu 



(fc) 



k=0 



where 9^ , and are the parameters of the models during the sub interval 
[t k ,t k+1 ] and A tk = t k+1 - t k . 

3 Fitting the Model 

To fit a stochastic volatility model to a given volatility surface, we first have to 
compute option prices, then to compute the implied volatility from these prices 
and finally to resolve the Optimization Problem. 

3.1 Optimization Problem 

We assume a set of option prices for different maturities and strikes as given. 
To calibrate a stochastic volatility model, we have to minimize an objective 
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function. We strive to minimize the mean squared error (MSE) of volatilities, 
i.e., 



minimize ^ (ffgiven - a fit) 2 , 

Options 

where <T g i V cn is a given vector of market Black-Scholes implied volatilities and am 
are the implied Black-Scholes volatility resulting from pricing the option with 
a stochastic volatility model. Instead of fitting the implied volatilities we could 
also fit option prices. Further it is possible to weight the differences between 
the given and the fitted prices. The choice of the weights may depend on the 
liquidity of a given option which can be assessed from the bid-ask spread. [Cont 



and Tankov 2004 Chapter 13] suggest to weight the option prices by squared 
Black-Scholes vegas evaluated at the implied volatilities of the market option 
prices. If we have a prior distribution, e.g. based on the posterior distribution 
from parameter estimation from historical data based on a Markov Chain Monte 
Carlo method, we may apply regularization techniques for which we refer to 



Cont and Tankov 2004 Chapter 13] 



3.2 Computation of Option prices 

For our application it is crucial to have an efficient method to price plain vanilla 
options. Therefore we use the Fast Fourier Transform (FFT) introduced for 
option pricing by Carr and Madan [1999] . The basic idea of this method is to 



develop an analytical expression for the Fourier transform of the option price 
and then to get the price back by Fourier inversion. Assuming no dividends and 
constant interest rates r, the initial value Co of a plain vanilla European call 
option is determined as 



Co = XoIT - Ke 



-rT 



n 2 



where 



IT = 



n 



2 — 



1 




/ Re 


7T 


'o 


1 


i r +oc 


2 + 


n Jo 



Re 



uX < K ) (f) T (u - i) 

iu4>t{—i) 
e- luln< - K )(f> T (u) 



du , 
du , 



with 4>t(u) being the characteristic function of \h(Xt)- Corresponding plain 
vanilla European put option prices can be derived via the Put-Call-parity equa- 
tion. Since we cannot evaluate the integral using the FFT, we rather use the 
expression 



r+oo 

C T (k) = / . 

Jk 



„-rT 



(e s - e k )q T (s)ds , 



where k denotes the log of the strike K and qr the risk-neutral density of the log 
price of the underlying. Since the call pricing function is not square integrable, 
we consider the modified call price cx(k) defined by 

&r(fc) = exp(afc)CV(fc) , 
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for a suitable a > 0. We consider now the Fourier transform of cy(fc) given by 



a 2 + a — v 2 + i(2a + l)v 
Then we can use the inverse transform to obtain Cr(fc) 

Qr (i) = 55t^re-*M»)*. 



» JO 

and the FFT to compute an approximation of CV(/t) via 

7T * — ' 

where 77 is the distance of the points of the integration grid and Vj = r)(j — 1). 



For more details we refer to Carr and Madan 1999 . This method is very fast 



since efficient implementations of the FFT exist. Moreover, this approach allows 
to compute the option prices for several strikes at once. Note that this approach 
only depends on the specific models via the characteristic function </>. 

3.3 Hybrid Optimization Algorithm 

We suggest a hybrid optimization algorithm which consists of a genetic algorithm 
with the purpose of locating a minimum within the specified bounds for the 
parameters and a subsequent direct search method to refine the solution. For 
a detailed reference concerning genetic algorithms we refer to Ba ck, Fogel, and| 



Michalewicz 20001. 



3.3.1 Genetic Algorithm (GA) 

Outline of the Algorithm The following outline summarizes how the genetic 
algorithm works: 

(i) The algorithm begins by creating randomly an initial population. The 
population is chosen such that possible constraints (e.g. upper or lower 
bounds on the variables) are satisfied. It is also possible to provide an ini- 
tial population or to specify only some members of the initial population. 

(ii) The algorithm evaluates the fitness of each individual in the population. 

(iii) The algorithm then creates iteratively a sequence of populations. At each 
step, the algorithm creates the next population based on the current gen- 
eration. Therefore the following steps are performed: 

(a) The algorithm evaluates the fitness of each individual in the popula- 
tion. The fitness is a measure on how well the optimization problem 
is solved. 

(b) A certain number of individuals in the current population (called the 
Elite count) that have the highest fitness values are chosen as elite 
children and are passed on to the next population. 



7 



(c) A certain number of children of the next generation are produced by 
combining a pair of parents (crossover) . The number is specified as 
fraction of the current population (Crossover fraction) . 

(d) The remaining children are produced by making random changes to 
a parent (Mutation). 

(e) The children of the current population form the next generation. 

(iv) The algorithm runs until one of the stopping criteria is met (e.g. maxi- 
mum number of (stall) generations/time reached, average change in fitness 
function less than a certain value, certain fitness value reached) 

For the Heston and Bates model we sample an initial population uniformly 
distributed between lower and upper bounds on the variables such that they 
satisfy the Feller condition. If parameter estimates from historical data are 
available these parameters may be included in the initial population. If the 
parameter estimates are based on a Markov Chain Monte Carlo method such 
that we get a posterior distribution for each parameter we may sample the initial 
population from the posterior distributions and the bounds on the parameters 
may be derived from the posterior distribution as well, e.g. as a quantile of the 
distribution. 

We use a crossover method that creates children as a linear combination of 
two parents, i.e., that lie on the line containing the two parents. We may choose 
the child randomly on the line between the two parents or alternatively a small 
distance away from the parent with the better fitness value in the direction 
away from the parent with the worse fitness value. Note that for the Heston 
model if parenti and parent2 with parameters (Oi, n i} pi,o-i,v ,i) satisfy the Feller 
condition 2/^0, — of > for i = 1,2, then the Feller condition is also satisfied 
for a parameter set on the line between the two parents since 

2(A«i + (1 - A)k 2 )(A0i + (1 - X)0 2 ) - (A<ti + (1 - A)cr 2 ) 2 > > 

for A e (0,1). 

For mutation we randomly generate directions that are adaptive with respect 
to the last successful or unsuccessful generation. A step length is chosen along 
each direction so that upper and lower bounds as well as the Feller constraint 
are satisfied. 

3.3.2 Direct Search Method 

After applying the GA with the purpose to locate a minimum in the specified 
bounds we use a pattern search algorithm to refine the result. We use complete 
polling to avoid local minima and only poll if the Feller constraint is satisfied. 

3.4 Fitting for Piecewise Constant Parameters 

For piecewise constant parameters we apply a bootstrapping method, i.e., after 
calibrating to one maturity we calibrate for the subsequent maturity, starting 
with the first maturity. We use as initial population for calibration the fitted 
parameters of the previous maturities. 



8 



Figure 1: given volatility surface 



Figure 2: fitted volatility surface 



If the calibration has to be done frequently, e.g. daily, we can achieve a fast 
recalibration by supplying additionally the parameters from the previous run as 
initial population for the new run (for the respective maturity). 

3.5 Numerical Example 

In this section we provide numerical results for fitting the Heston model. We con- 
sider an FX volatility surface for the Yen/Dollar exchange rate from 21/05/2008 
with values <7 g i ve n for a given set of delta values A g i vcn . Since we need the volatil- 
ity surface in terms of strikes, we transform A g i ven into K g i vcn using 



K ■ 

1 kgivcn 



s 



CXp^- 1 (Agivcn)^ - ( r + s_)t) 



where Af denotes the standard normal cumulative distribution function. More- 
over, for at the money options, we have to substitute the strike with the cor- 
responding forward price. To obtain crgt, we compute the option prices and 
retrieve the corresponding Black-Scholes model implied volatilities cr^ t via a 
standard bisection method. For the computation of option prices we use the 



FFT as described in Section 3.2 We use a = 0.75 since Schoutens et al. 2004 
show in an empirical study that with this value prices are well replicated for 
many model parameters. The parameters N — 2 16 and i] = 0.015 turned out to 
be a reasonable trade-off between accuracy and speed. 

Figure [T] shows the given volatility surface quoted in terms of Deltas and 
Figure[2]is the surface produced by our fitted model. Figure|3]gives the difference 
between the quoted surface and the fitted surface. Figure [4] illustrates how 
the parameters evolve over time. Due due the fact that the calibration is an 
ill posed problem, i.e. the solution might not be unique and the parameters 
don't depend continuously on the data, we see considerable jumps of the time 
dependent parameters. Figure [5] fl2] show the given smile and the fitted smile 
for the various maturities. 



4 Conclusion and Outlook 



We get a good fit of the given volatility surface for stochastic volatility mod- 
els using piecewise constant parameters using a hybrid optimization algorithm 




Time to maturity 



0.08 
0.06 
0.04 
0.02 



3A66M9M1Y 2Y 

Time to maturity 



3R66M 9M 1 Y 2Y 3Y 3H66M 9M 1 Y 2Y 

Time to maturity Time to maturity 

Figure 4: (time dependent) parameters 
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Figure 5: fit for maturity 3M 



Figure 6: fit for maturity 4M 




and can use this fit to price exotic derivatives. We also note that the piece- 
wise constant parameters may jump considerably. This is due to the fact the 
the calibration problem per se is ill posed giving evidence to the necessity of 
regularization techniques. If the posterior distribution from parameter estima- 
tion from historical data is available, we may apply regularization techniques. 
Further, the bounds for the parameters for the optimization algorithm could 
be derived from (quantiles of) the posterior distribution and also the initial 
population could be sampled from the posterior distribution. 
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